Python/NEURON code for simulating biophysically realistic thalamocortical dynamics during sleep

Understanding the function of sleep and its associated neural rhythms is an important goal in neuroscience. While many theoretical models of neural dynamics during sleep exist, few include the effects of neuromodulators on sleep oscillations and describe transitions between sleep and wake states or different sleep stages. Here, we started with a C++-based thalamocortical network model that describes characteristic thalamic and cortical oscillations specific to sleep. This model, which includes a biophysically realistic description of intrinsic and synaptic channels, allows for testing the effects of different neuromodulators, intrinsic cell properties, and synaptic connectivity on neural dynamics during sleep. We present a complete reimplementation of this previously-published sleep model in the standardized NEURON/Python framework, making it more accessible to the wider scientific community.


Introduction
The phenomenon of sleep is mysterious: nearly all animals exhibit some form of sleep, yet we cannot say exactly why this must be so.Many clues have surfaced over the past few decades, with neural activity during mammalian sleep exhibiting a rich repertoire of rhythms associated with memory consolidation [1,2], toxin clearance [3], and regulation of synaptic strength [4].Different stages of sleep feature characteristic neural activity patterns, with sleep spindles prominent during non-REM Stage 2 (associated with memory consolidation [5]), slow waves dominating non-REM Stage 3 (associated with memory consolidation [6] and clearance of metabolic waste [7]), and wake-like alpha activity characteristic of REM sleep.Such neural activity is heavily influenced by the concentration of various neuromodulators (such as acetylcholine, histamine, norepinephrine, etc.) that modulate the activity of neurons at the cellular and network levels.
Disentangling the effects of various neuromodulators on neural activity, and how this influences sleep activity and transitions between sleep stages, is an active area of exploration in neuroscience.Experimental investigation of these questions is ongoing [8-10], but it is often difficult or impossible to determine levels of neuromodulators across the brain with current technology.Computationally modeling the effects of neuromodulators on neural activity during sleep can provide useful models to guide experimental investigation.One such computational model that examines the effect of neuromodulation on sleep activity was proposed by Krishnan et al. [11].It simulates the effects of varying concentrations in acetylcholine, histamine, and GABA on neural activity during wake, N2, N3, and REM sleep.The model was based on extensive previous modeling and experimental works [12][13][14][15][16][17][18][19][20][21][22][23][24], and several subsequent studies have utilized this model's approach to simulate neuromodulation in their work [25][26][27][28][29][30].It is also one of the few biophysically realistic, detailed individual neuron-based (as opposed to mean field or firing rate-based) models of neural dynamics during sleep.We have re-implemented it in a more accessible framework that allows other researchers to easily modify the model to explore their own lines of enquiry.

Software
The model involves a simulated thalamocortical loop, with 500 cortical pyramidal (PY) neurons, 100 cortical inhibitory (INH) neurons, 100 thalamocortical (TC) relay neurons, and 100 thalamic reticular nucleus (RE) neurons.Krishnan et al. [11] used this model to show how plausible changes in neuromodulation could give rise to the neural rhythms characteristic of different stages of sleep, and they used their results to suggest that differences in neural rhythms between human and animal recordings were primarily due to differences in acetylcholine concentration.The model's code was originally written in C++, which resulted in fast simulation times, but it is hardly accessible for newcomers and results in slow development of new features.To ease this situation we present here the model from [11], ported to the widely used Python-based NEURON simulation environment [31].Fig. 1 depicts the main simulation results from [11], but regenerated using NEURON.NEURON is one of the leading simulation environments in the field of computational neuroscience, allowing users to develop models at the level of biophysical detail using Python, without having to worry about numerical implementation.NEURON auto-generates C code to take care of such numerical details, so that simulations are fast without requiring inordinate development time.NEURON also includes libraries for parallel computation, and our code is fully parallelized [32].The 360,000 ms simulation shown in Fig. 1 took approximately 15 min to run with 128 cores on the Neuroscience Gateway Computing (NSG) cluster [33].Commands for submitting this simulation to NSG are included in the repository.Further, the model is also available at the Cybershuttle site (https://neuroscience.cybershuttle.org)[34], which allows for tracking simulations and parameter sweeps.
The code is organized around four Python files: config.pydefines all simulation parameters, cell_classes.pydefines classes for the four kinds of neurons (PY, INH, TC, and RE), network_class.pydefines a class prescribing network connectivity, and main.pyinitializes and runs the simulation.The most important ancillary files are the so-called "mod files," which consist of high-level NMODL code specifying various mathematical details of the model (such as equations describing the biophysics of various kinds of synapses).mod files are translated into C code by NEURON before running a simulation (using the command nrnivmodl), which is how NEURON simulations can run so efficiently without the need to explicitly develop code in a low-level language like C or C++.In addition, two legacy files are included (interpxyz.hoc and calcrxc_a.hoc)for simulations in which the cortical LFP is calculated from biophysical first principles.These compute the spatial locations of all neuronal compartments and the simulated recording electrode, as well as the transfer resistance between the recording electrode and each compartment.The detailed instructions for installation and running are described in the README.mdfile of the software package.

Impact
Our NEURON simulation code serves as a foundation for other researchers interested in developing computational models of thalamocortical activity during wake and sleep.Additional capabilities such as simulating stochastically varying neuromodulator concentrations, memory consolidation, spindle-slow wave coupling, etc., are all much more readily developed in NEURON than in lower-level languages such as C++.One feature that we have added as a new feature in this software package is the calculation of local field potential (LFP) from biophysical first principles, as opposed to simply summing intracellular voltages.While summing intracellular voltages is a useful first approximation to the biophysical LFP, this method does not reflect the true source of extracellular potentials.Our NEURON implementation uses the method outlined in [35] to determine the contributions of all transmembrane currents to the net potential at the recording electrode.Given the widespread interest in how different sleep rhythms emerge and interact in different stages of sleep [36][37][38], the more realistic simulation of these rhythms offered by our model may be of significant interest to the broader scientific community.Other researchers may also be interested in using our "mod" files in their own NEURON models.For example, the mod file describing synaptic interactions between pyramidal cells (AMPA_D2.mod) is essential to the generation of cortical slow waves in our model, and it could be used in other cortical models for the same purpose.
More broadly, previous work has suggested that simulated sleep can improve deep learning paradigms.One flaw in modern artificial neural network performance is the phenomenon of "catastrophic forgetting," in which the learning of new memories invariably degrades previously-learned memories.The human brain does not suffer from this problem to nearly the same extent, and it is thought that sleep is a major reason why.This has inspired efforts to simulate sleep in artificial neural networks, in order to mitigate catastrophic forgetting [29,39].Our Python/NEURON implementation of this model of sleep will enable more nimble development of software to address the problem of catastrophic forgetting.
In the future, we plan to develop this model to simulate the effects of various neuromodulators on sleep rhythms [40][41][42][43].Future development could include upgrading the code to be compatible with coreNEURON [44], which would allow it to run on GPU architectures.

Limitations
From a neuroscience point of view, this thalamocortical model is limited by the fact that it includes the effects of only a select group of neuromodulators (omitting norepinephrine, orexin, and serotonin, among others).However, incorporating these neuromodulators into the model would be relatively straightforward using NEURON.This network model also uses relatively simple one-or two-compartment models for individual neurons.The original C++ code further simplified the two-compartment cortical cells by omitting capacitance from the axosomatic compartment, which helped to speed the simulation numerically [12,45].In porting the model to NEURON, we included capacitance in the axosomatic compartment, which required that associated conductances be downscaled by a factor of 3/4 in order to give similar qualitative results to the C++ code.The NEURON code is therefore not an exact translation of the C++ code, but comparing Fig. 1 of this paper to Fig. 2 in [11] demonstrates the overall fidelity of the NEURON code to the original.

Conclusion
The thalamocortical model by Krishnan et al. [11] is one of the few biophysically realistic network models of different sleep stages in the computational neuroscience literature.
Here we have re-written the code using the Python-based NEURON framework, enabling other researchers to efficiently adapt the model to their own research questions.Anyone interested in computational modeling of sleep rhythms, in particular their dependence on neuromodulatory tone and their possible effects on memory consolidation, will find this software useful.

Fig. 1 .
Fig. 1.Simulated sleep cycles In a thalamocortical model.(A) Raster plot of network activity for wake, N2, N3, REM, and N2 sleep stages.Different sleep stages were generated my modulating acetylcholine, histamine, and GABA.(B) Spectrogram of simulated cortical local field potential.Note the spindle rhythms in N2 and the slow waves in N3. (C) Simulated cortical local field potential, computed by summing intracellular voltage traces.